;;..................................................................
;; zonal aerosol water content 
;;..................................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    

;;..................................................................
;; interp_pres
;;..................................................................
  undef("interp_pres_cam")

  function interp_pres_cam(va,PS)

  begin

  L_debug = "1" 

  nlat = dimsizes(va(0,0,:,0))    
  nlon = dimsizes(va(0,0,0,:))   
  nlev = dimsizes(va(0,:,0,0))   
  nt = 1 

  I_lev_14 = 0
  I_lev_10 = 1

  if(I_lev_14.eq.0) then
     lev = new(14,"float")
     lev = (/10,20,50,70,100,150,200,300,400,500,600,700,850,1000/)
  end if

  if(I_lev_10.eq.0) then
     lev = new(10,"float")
     lev = (/100,200,300,400,500,600,700,800,900,1000/)
  end if


  if(I_lev_14.eq.0) then
     vp = new( (/ nt, 14, nlat, nlon /), float )
  end if

  if(I_lev_10.eq.0) then
     vp = new( (/ nt, 10, nlat, nlon /), float )
  end if

  vp!0 = "time" 
  vp&time = va&time 
  vp!1 = "lev" 
  vp&lev = lev
  vp&lev@units = "hPa" 
  vp!2 = "lat" 
  vp&lat = va&lat
  vp!3 = "lon" 
  vp&lon = va&lon

  pold  = new( (/nlev,nlat,nlon/), float ) ;hpa

;;..................................................................
;; load vertical coordinate
;;..................................................................
;;load "mkdata_mm_vpara.ncl"

 hyam_cam = (/ 0.00364346569404006, 0.00759481964632869, 0.0143566322512925, \  
    0.0246122200042009, 0.0382682997733355, 0.0545954797416925,       \  
    0.0720124505460262, 0.0878212302923203, 0.103317126631737,        \  
    0.121547240763903, 0.142994038760662, 0.168225079774857,          \  
    0.178230673074722, 0.170324325561523, 0.161022908985615,          \  
    0.150080285966396, 0.137206859886646, 0.122061938047409,          \  
    0.104244712740183, 0.0849791541695595, 0.0665016956627369,        \  
    0.0501967892050743, 0.037188658490777, 0.028431948274374,         \ 
    0.0222089774906635, 0.016407382208854, 0.0110745579004288,        \     
    0.00625495356507599, 0.00198940909467638, 0                       /) 

 hybm_cam = (/ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0196774136275053,  \  
               0.062504293397069, 0.112887907773256, 0.172161616384983,    \  
               0.241894043982029, 0.323930636048317, 0.420442461967468,    \  
               0.524799540638924, 0.624887734651566, 0.713207691907883,    \  
               0.783669710159302, 0.831102818250656, 0.864811271429062,    \  
               0.896237164735794, 0.92512384057045, 0.951230525970459,     \  
               0.974335998296738, 0.992556095123291                        /) 

;;..................................................................
;; pressure (hPa) 
;;..................................................................
  do it = 0, nt - 1
    do k = 0, nlev - 1
      pold(k,:,:) = ( PS(:,:)/100.) * (hybm_cam(k)) + 1013.25 * (hyam_cam(k))
      ;;pold(k,:,:) = ( PS(:,:)/100.) * (hybm_echam(k) ) + ( hyam_echam(k) ) / 100.
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do
  end do 


;;..................................................................
;; veritcal interpolation
;;..................................................................
  L_vinth2p = True

  if (L_vinth2p) then
     ;;fake = hyam_echam
     ;;tmpp = vinth2p(varx(it,:,:,:),fake,hybm_echam,lev*100.,PS(:,:)*100.,1,1.,1,False)
     ;;fake = hyam_echam / 101325.
     ;;tmpp = vinth2p(va(0,:,:,:),fake,hybm_echam,lev,PS(:,:),1,1013.25,1,False)

     fake = hyam_cam
     tmpp = vinth2p(va(0,:,:,:),fake,hybm_cam,lev,PS(:,:),1,1000.,1,False)


     ;;print("max,min va: "+max(va(0,:,:,:)) +" "+min(va(0,:,:,:)))
     ;;print("max,min tmpp: "+max(tmpp) +" "+min(tmpp))

     vp(0,:,:,:) = (/ tmpp(:,:,:) /)

    ;printVarSummary(vp)
    ;printVarSummary(tmpp)

  else

    do   j = 0, nlat - 1
      do i = 0, nlon - 1
        vp(0,:,j,i) = linint1( pold(:,j,i), va(0,:,j,i), False, lev, 0 )
      end do
    end do

  end if

  do i=0,5
     vp(0,i,:,:) = mask(vp(0,i,:,:),(lev(i).gt.PS(:,:)/100.),False)
  end do

  return(vp)

end

;;..................................................................
;; interp_pres
;;..................................................................
  undef("interp_pres_echam")

  function interp_pres_echam(va,PS)

  begin

  L_debug = "1" 

  nlat = dimsizes(va(0,0,:,0))    
  nlon = dimsizes(va(0,0,0,:))   
  nlev = dimsizes(va(0,:,0,0))   
  nt = 1 

  I_lev_14 = 0
  I_lev_10 = 1

  if(I_lev_14.eq.0) then
     lev = new(14,"float")
     lev = (/10,20,50,70,100,150,200,300,400,500,600,700,850,1000/)
  end if

  if(I_lev_10.eq.0) then
     lev = new(10,"float")
     lev = (/100,200,300,400,500,600,700,800,900,1000/)
  end if


  if(I_lev_14.eq.0) then
     vp = new( (/ nt, 14, nlat, nlon /), float )
  end if

  if(I_lev_10.eq.0) then
     vp = new( (/ nt, 10, nlat, nlon /), float )
  end if

  vp!0 = "time" 
  vp&time = va&time 
  vp!1 = "lev" 
  vp&lev = lev
  vp&lev@units = "hPa" 
  vp!2 = "lat" 
  vp&lat = va&lat
  vp!3 = "lon" 
  vp&lon = va&lon

  pold  = new( (/nlev,nlat,nlon/), float ) ;hpa

;;..................................................................
;; load vertical coordinate
;;..................................................................
  load "mkdata_mm_vpara.ncl"

;;..................................................................
;; pressure (hPa) 
;;..................................................................
  do it = 0, nt - 1
    do k = 0, nlev - 1
      pold(k,:,:) = ( PS(:,:)/100.) * (hybm_echam(k) ) + ( hyam_echam(k) ) / 100.
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do
  end do 


;;..................................................................
;; veritcal interpolation
;;..................................................................
  L_vinth2p = True

  if (L_vinth2p) then
     ;;fake = hyam_echam
     ;;tmpp = vinth2p(varx(it,:,:,:),fake,hybm_echam,lev*100.,PS(:,:)*100.,1,1.,1,False)
     fake = hyam_echam / 101325.
     tmpp = vinth2p(va(0,:,:,:),fake,hybm_echam,lev,PS(:,:),1,1013.25,1,False)

     ;;print("max,min va: "+max(va(0,:,:,:)) +" "+min(va(0,:,:,:)))
     ;;print("max,min tmpp: "+max(tmpp) +" "+min(tmpp))

     vp(0,:,:,:) = (/ tmpp(:,:,:) /)

    ;printVarSummary(vp)
    ;printVarSummary(tmpp)

  else

    do   j = 0, nlat - 1
      do i = 0, nlon - 1
        vp(0,:,j,i) = linint1( pold(:,j,i), va(0,:,j,i), False, lev, 0 )
      end do
    end do

  end if

  do i=0,5
     vp(0,i,:,:) = mask(vp(0,i,:,:),(lev(i).gt.PS(:,:)/100.),False)
  end do

  return(vp)

end



;;..................................................................
;; main program 
;;..................................................................

begin 

;;..................................................................
;; plot configuration 
;;..................................................................
  np = 1

  fontscale = 1.2 

;;..................................................................
;; read data 
;;..................................................................

  fnbase = "X0006A_0003-0007_ANN" 

  fna = "wsubmm_T63L31K01_2010_ANN.nc" 
  fns = "metmm_T63L31K01_2010_ANN.nc" 
  fna = "wsubmm_X0006A_0003-0007_ANN.nc" 
  fns = "PS_X0006A_0003-0007_ANN.nc"

  fla = addfile(fna,"r")
  fls = addfile(fns,"r")

  va = fla->WSUBO ;;W_TURB
  aps= fls->PS    ;;aps

  PS = dim_avg_Wrap(aps(lat|:,lon|:,time|:))

  nlat = dimsizes(va(0,0,:,0))
  nlon = dimsizes(va(0,0,0,:))
  nlev = dimsizes(va(0,:,0,0))   
  nt = 1

;;..................................................................
;; interpolation to pressure leves 
;;..................................................................

  model = "echam" 
  model = "cam" 

  if(model.eq."echam") then 
     vpa = interp_pres_echam(va,PS) 
  else 
     vpa = interp_pres_cam(va,PS) 
  end if 

;;..................................................................
;; zonal average 
;;..................................................................
  vza = dim_avg_Wrap(vpa(lev|:,lat|:,time|:,lon|:))  

;; lev = (/10,20,50,70,100,150,200,300,400,500,600,700,850,1000/)

  i200 = 6 
  i300 = 7 
  i500 = 9 
  i700 = 11
  i850 = 12

;;..................................................................
;; 200 hPa 
;;..................................................................
  v200 = vpa(:,i200,:,:) 

;;..................................................................
;; 300 hPa 
;;..................................................................
  v300 = vpa(:,i300,:,:) 

;;..................................................................
;; 500 hPa 
;;..................................................................
  v500 = vpa(:,i500,:,:) 

;;..................................................................
;; 700 hPa 
;;..................................................................
  v700 = vpa(:,i700,:,:) 

;;..................................................................
;; 850 hPa 
;;..................................................................
  v850 = vpa(:,i850,:,:) 

;;..................................................................
;; output data 
;;..................................................................

  fno = fnbase + "_pres.nc"  

  system("rm "+fno) 

  flo = addfile(fno,"c") 

  flo->v200=v200 
  flo->v300=v300 
  flo->v500=v500 
  flo->v700=v700 
  flo->v850=v850 
 

  fno = fnbase + "_zonal.nc"  

  system("rm "+fno) 

  flo = addfile(fno,"c") 

  flo->vz=vza

;;..................................................................
;; colormap 
;;..................................................................
  colormap = "WhiteBlueGreenYellowRed" 

;;..................................................................
;; open wks 
;;..................................................................
  wks = gsn_open_wks("eps","out")            ; open a ps file
  gsn_define_colormap(wks,colormap)      ; choose colormap

;;..................................................................
;; plot configure
;;..................................................................
  plot = new(np,graphic)                       ; create a plot array
  
;;..................................................................
;; plot defination 
;;..................................................................
  res                     = True
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnInfoLabelOn       = False           ; turn off cn info label
  res@cnFillOn            = True            ; turn on color
  res@cnLinesOn            = False              ; turn off contour lines
  res@gsnSpreadColors     = True            ; spread out color table
    
  res@gsnSpreadColorStart = 4
  ;;res@gsnSpreadColorEnd   = -3 
  
  res@lbOrientation    =   "vertical"          ; vertical label bar
  res@lbLabelFontHeightF = 0.021 * fontscale 
  res@lbTitleOn        = False
  res@lbTitleString    = "" 
  res@lbTitleFontHeightF= .015                 ; make title smaller
  ;;res@lbLabelAngleF    =    45
  ;;res@pmLabelBarParallelPosF   = 0.00 
  ;;res@pmLabelBarOrthogonalPosF = 0.12           ; move whole thing down

  res@vpWidthF             = 0.8        ; width of contour plots

  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels            = (/0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,1.0,1.2,1.5/) 

;;  res@cnLevelSelectionMode =  "ManualLevels"   
;;  res@cnMinLevelValF       = -10.
;;  res@cnMaxLevelValF       =  45.
;;  res@cnLevelSpacingF      =   5. 

;;..................................................................
;; Font 
;;..................................................................
  res@gsnStringFontHeightF = 0.025 * fontscale    
  res@tmXBLabelFontHeightF = 0.022 * fontscale  
  res@tmYLLabelFontHeightF = 0.022 * fontscale  
  res@tiMainFontHeightF    = 0.025 * fontscale    
  res@tiXAxisFontHeightF   = 0.025 * fontscale  
  res@tiYAxisFontHeightF   = 0.025 * fontscale  

;;..................................................................
;; individual plot 
;;..................................................................

  res@gsnRightString = "kg kg~S~-1~N~" 
  res@gsnRightString = "m s~S~-1~N~" 

  ip = 0 

  ks = 4 ;;6

  res@gsnLeftString = "Charasteristic updraft velocity" 
  res@tiXAxisString = "" 
  res@tiYAxisString = "Pressure (hPa)" 
  plot(ip) = gsn_csm_pres_hgt(wks,vza(ks:,:,0),res) 
  ip = ip + 1 

;;;  res@tiYAxisString = "" 
;;;  res@gsnLeftString = "b) HAM2 default (~F33~k~F~-Koehler 95%)" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzb(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;  res@tiXAxisString = "Latitude"  
;;;  res@tiYAxisString = "Pressure (hPa)" 
;;;  res@gsnLeftString = "c) HAM2 with ZSR 100%" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzc(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;  res@tiYAxisString = "" 
;;;  res@gsnLeftString = "d) HAM2 with ZSR 95%" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzd(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;;;..................................................................
;;;;; global mean aod, wataod, burden 
;;;;;
;;;;;   cases        ttaod     wataod     burden 
;;;;; 
;;;;;   HAM1         0.140      0.105       75.2
;;;;;   HAM2         0.135      0.094       50.9
;;;;;   ZSR 100%     0.163      0.123       99.9
;;;;;   ZSR  95%     0.138      0.097       56.0 
;;;;;..................................................................
;;;
;;;  ttaod = (/0.140,0.135,0.163,0.138/)
;;;  awaod = (/0.105,0.094,0.123,0.097/)
;;;  awpct = (/74.8,69.8,75.1,69.7/)
;;;  burden= (/75.2,50.9,99.9,56.0/)
;;;
;;;;;..................................................................
;;;;; add text : global mean aod, wataod, burden 
;;;;;..................................................................
;;;  txres               = True
;;;  txres@txFontHeightF = 0.022 * fontscale 
;;;  txres@Position      = "Left"
;;;  txres@txJust        = "CenterLeft"
;;;
;;;  texta = new(np,graphic)
;;;  textb = new(np,graphic)
;;;  textc = new(np,graphic)
;;;  
;;;  do ip = 0,np-1
;;;
;;;     ;;ca = "Total AOD        " 
;;;     cb = "Water AOD      " 
;;;     cc = "Water Burden   " 
;;;
;;;     ;;printa = ca + sprintf("%3.3f",ttaod(ip)) 
;;;     printb = cb + sprintf("%3.3f",awaod(ip)) + " ("+sprintf("%3.1f",awpct(ip))+"%)"  
;;;     printc = cc + sprintf("%3.1f",burden(ip))+" Tg" 
;;;
;;;     ;;texta(ip) = gsn_add_text(wks,plot(ip),printa,-60.,250.,txres)
;;;     textb(ip) = gsn_add_text(wks,plot(ip),printb,-60.,340.,txres)
;;;     textc(ip) = gsn_add_text(wks,plot(ip),printc,-60.,290.,txres)
;;;
;;;  end do

;;..................................................................
;; create panel
;;..................................................................
  resP                  = True                   ; modify the panel plot
  resP@gsnPanelLabelBar = False 
  ;;resP@gsnPanelBottom   = 0.05
  resP@txString   = ""
  resP@gsnPanelYWhiteSpacePercent = 5
  resP@gsnPanelXWhiteSpacePercent = 0
  gsn_panel(wks,plot,(/2,2/),resP)               ; now draw as one plot


end 









